Circle

Go back to the About page.

Let us set some global options for all code chunks in this document.

knitr::opts_chunk$set(
  message = FALSE,    # Disable messages printed by R code chunks
  warning = FALSE,    # Disable warnings printed by R code chunks
  echo = TRUE,        # Show R code within code chunks in output
  include = TRUE,     # Include both R code and its results in output
  eval = TRUE,       # Evaluate R code chunks
  cache = FALSE,       # Enable caching of R code chunks for faster rendering
  fig.align = "center",
  #out.width = "100%",
  retina = 2,
  error = TRUE,
  collapse = FALSE
)
rm(list = ls())
set.seed(1982)

1 Import libraries

# Install R-INLA package
# install.packages("INLA",repos = c(getOption("repos"),INLA ="https://inla.r-inla-download.org/R/testing"), dep = TRUE)
# Update R-INLA package
# inla.upgrade(testing = TRUE)
# Install inlabru package
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# Install rSPDE package
# remotes::install_github("davidbolin/rspde", ref = "devel")
# Install MetricGraph package
# remotes::install_github("davidbolin/metricgraph", ref = "devel")

library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(Matrix)

library(dplyr)
library(plotly)
library(scales)
library(patchwork)
library(tidyr)
library(ggplot2)
library(reshape2)
library(sf)

library(here)
library(rmarkdown)
library(knitr)
library(grateful) # Cite all loaded packages

library(latex2exp)
library(plotrix)

rm(list = ls()) # Clear the workspace
set.seed(1982) # Set seed for reproducibility

2 Circle graph

Follow this link for more details.

# Matern covariance function
matern.covariance <- function(h, kappa, nu, sigma) {
  if (nu == 1 / 2) {
    C <- sigma^2 * exp(-kappa * abs(h))
  } else {
    C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
      ((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
  }
  C[h == 0] <- sigma^2
  return(as.matrix(C))
}

# Folded.matern.covariance.1d
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma, L = 1, N = 10, boundary = c("neumann",
                                                                                               "dirichlet", "periodic")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
    stop("The possible boundary conditions are 'neumann',
    'dirichlet' or 'periodic'!")
  }
  addi = t(outer(x, x, "+"))
  diff = t(outer(x, x, "-"))
  s1 <- sapply(-N:N, function(j) { # s1 is a matrix of size length(h)x(2N+1)
    diff + 2 * j * L
  })
  s2 <- sapply(-N:N, function(j) {
    addi + 2 * j * L
  })
  if (boundary == "neumann") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) +
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else if (boundary == "dirichlet") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) -
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else {
    C <- rowSums(matern.covariance(h = s1,
                                   kappa = kappa, nu = nu, sigma = sigma))
  }
  return(matrix(C, nrow = length(x)))
}

# Function to get the true covariance matrix
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
  h = c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[,2])
  true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
  return(true_cov_mat)
}

# Define the graph
r <- 1/(pi)
theta <- seq(from = -pi, to = pi, length.out = 100)
edge <- cbind(1+r+r*cos(theta), r*sin(theta))
edges <- list(edge)
graph <- metric_graph$new(edges = edges)

# parameters
h.ok <- 2^-10
type <- "covariance"
type_rational_approximation = "brasil"
rho <- 0.5
#m = 4
sigma <- 1
N.folded <- 10
boundary <- "periodic" # Do not change this

# Mesh sizes
h_aux <- seq(5.5, 4, by = -1/2)
h_vector <- 2^-h_aux
h_label <- paste0("2^-", h_aux, "")
h_label_latex <- sprintf("$2^{-%f}$", h_aux)

# Beta values
beta_aux <- c(4, 6, 6.5, 7, 7.5)
beta_vector <- beta_aux/8
beta_label <- paste0(beta_aux, "/8")
theoretical_rate <- pmin(4*beta_vector-1/2,2)

graph.ok <- graph$clone()
# Build graph with overkill mesh
graph.ok$build_mesh(h = h.ok)
graph.ok$compute_fem()

# Get the overkill mesh locations
loc.ok <- graph.ok$mesh$VtE # or graph.ok$get_mesh_locations()

# Initialize the list of graphs and the list of projection matrices
graphs <- list()
A <- list()
for(i in 1:length(h_vector)){
  graphs[[i]] <- graph$clone()
  graphs[[i]]$build_mesh(h = h_vector[i])
  A[[i]] <- graphs[[i]]$fem_basis(loc.ok)
}

cov.error.folded <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))

m_values <- c()
for (j in 1:length(beta_vector)) {
  beta <- beta_vector[j]
  alpha <- 2*beta
  fract2beta <- alpha - floor(alpha)
  nu <- alpha - 0.5
  kappa <- sqrt(8*nu)/rho
  tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1

  Sigma.folded <- gets_true_cov_mat(graph = graph.ok, 
                               kappa = kappa,
                               nu = nu,
                               sigma = sigma,
                               N = N.folded,
                               boundary = boundary)

  for (i in 1:length(h_vector)) {
    h <- h_vector[i]
    m <- min(20, ceiling((min(4*beta - 1/2,2) + 1/2)^2*log(h)^2/(4*pi^2*fract2beta)))
    m_values <- c(m_values, m)

    Sigma <- matern.operators(alpha = alpha,  
                             kappa = kappa,
                             tau = tau,
                             m = m,
                             graph = graphs[[i]],
                             type = type,
                             type_rational_approximation = type_rational_approximation)$covariance_mesh()

    Sigma.approx <- A[[i]]%*%Sigma%*%t(A[[i]])

    cov.error.folded[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.folded - Sigma.approx)^2%*%graph.ok$mesh$weights))
  }
}


print(m_values)

##  [1] 20 20 20 20  5  4  4  3  4  4  3  2  4  3  3  2  3  3  2  2

slope_circle <- numeric(length(beta_vector))

for (u in 1:length(beta_vector)) {
  slope_circle[u] <- coef(lm(log(cov.error.folded[,u]) ~ log(h_vector)))[2]
}

transposed_df <- data.frame(t(data.frame(beta = beta_label, theoretical_rate = theoretical_rate, slope = slope_circle)))
rownames(transposed_df) <- c("beta", "Theoretical rates", "Circle graph")
colnames(transposed_df) <- NULL
# Display the transposed data frame
transposed_df |> paged_table()

loglog_line_equation <- function(x1, y1, slope) {
  b <- log10(y1 / (x1 ^ slope))

  function(x) {
    (x ^ slope) * (10 ^ b)
  }
}

guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
for (j in 1:length(beta_vector)) {
  guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
  for(k in 1:length(h_vector)){
    point_x1 <- h_vector[k]
    point_y1 <- cov.error.folded[k, j]
    slope <- theoretical_rate[j]
    line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
    guiding_lines_aux[,k] <- line(h_vector)
  }
  guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}

# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))

# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
  geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
            aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})

df <- as.data.frame(cbind(h_vector, cov.error.folded))
colnames(df) <- c("h_vector", beta_label)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")

p2 <- ggplot() +
  geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  plot_lines +
  labs(title = "Circle graph",
    x = TeX("$\\hat{h}$"),
       y = "Covariance Error",
       color = expression(beta)) +
  scale_x_log10(breaks = h_vector, labels = TeX(h_label_latex)) +
  scale_y_log10() +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

p2